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Abstract 

We  study  the  discretization  behavior  of  classical  finite  element  and  NURBS  approxima¬ 
tions  on  problems  of  structural  vibrations  and  wave  propagation.  We  find  that,  on  the  basis 
of  equal  numbers  of  degrees-of-freedonr  and  bandwidth,  NURBS  have  superior  approxima¬ 
tion  properties.  In  fact,  we  observe  that  the  high  mode  behavior  of  classical  finite  elements 
is  divergent  with  the  order  of  approximation,  a  surprisingly  negative  result.  On  the  other 
hand,  NURBS  offer  almost  spectral  approximation  properties,  and  all  modes  converge  with 
increasing  order  of  approximation. 
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1  Introduction 


Isogeometric  Analysis  was  introduced  by  Hughes,  Cottrell  and  Bazilevs  [15]  in  an 
effort  to  improve  upon  shortcomings  of  finite  element  analysis  in  the  areas  of  ge¬ 
ometric  precision,  ease  of  mesh  refinement,  and  integration  with  Computer  Aided 
Design  (CAD).  Engineering  investigations  of  performance  in  vibration  calculations 
showed  very  good  behavior,  Cottrell  et  al.  [7],  as  did  mathematical  studies  of  con¬ 
vergence  under  mesh  refinement,  Bazilevs  et  al.  [2].  Nevertheless,  the  approxima- 
bility  of  isogeometric  analysis  compared  with  classical  finite  element  analysis  has 
not  been  thoroughly  investigated.  It  is  the  purpose  of  this  paper  to  initiate  such  a 
comparison. 

To  a  certain  degree,  one  might  say  that  isogeometric  analysis  subsumes  finite  ele¬ 
ment  analysis,  in  that  standard  element  basis  functions  can  be  generated  from  B- 
splines  and  NURBS  (non-uniform  rational  B-splines),  the  main  technologies  used 
thus  far  in  the  instantiation  of  isogeometric  analysis.  However,  isogeometric  anal¬ 
ysis  offers  more,  in  particular,  the  possibility  of  developing  smoother  basis  func¬ 
tions,  at  least  throughout  patches  (i.e.,  subdomains),  and  in  many  cases  globally. 
Isogeometric  analysis  emanates  from  constructs  used  in  design,  computer  graph¬ 
ics,  animation,  and  visualization,  and  it  is  often  the  case  that  smoothness  is  of  ut¬ 
most  importance.  For  example,  rendering  of  reflective  objects  requires  essentially 
C'2-continuity  in  order  not  to  exhibit  spurious  reflections.  In  recent  years,  a  num¬ 
ber  of  computer  graphics  techniques  have  been  developed  to  address  this  issue. 
It  is  interesting  to  note  that  difficulties  in  developing  smooth  basis  functions  (i.e., 
G'1 -continuity  and  higher)  in  the  early  years  of  finite  elements  usually  led  to  re¬ 
formulation  of  problems  so  that  Cy°-elements  could  be  utilized,  a  classical  illustra¬ 
tion  being  the  virtual  abandonment  of  Poisson-Kirchhoff  plate  theory,  which  leads 
to  fourth-order  biharmonic  problems,  in  favor  of  Reissner-Mindlin  theory,  which 
leads  to  second-order  differential  equations.  Times  have  changed  and  it  is  now  pos¬ 
sible  to  construct  complex  models  utilizing  functions  smoother  than  C°.  This  may 
open  a  door  to  simpler  formulations  of  problems  involving  higher-order  differen¬ 
tial  operators  (see,  e.g.,  Gomez  et  al.  [11]).  Another  possibility  is  that  smoother 
functions  might  produce  better  approximations  of  derivatives  than  C°  continuous 
finite  elements  in  second-order  problems.  For  example,  stresses  are  generally  the 
most  important  quantities  in  structural  analysis,  and  they  are  usually  smooth  al¬ 
most  everywhere.  Standard  finite  elements  require  smoothing  and  post-processing 
of  stresses.  This  might  be  avoided  through  the  use  of  smooth  basis  functions. 

In  this  paper  we  initiate  the  investigation  of  smooth  basis  functions  generated  by 
isogeometric  analysis  and  compare  them  with  standard  C°  finite  elements.  The 
problems  used  as  a  basis  of  comparison  emanate  from  structural  dynamics  and 
wave  propagation,  in  particular  the  eigenvalue  problem  of  free  vibration,  and  the 
Helmholtz  equation  of  time  harmonic  wave  propagation  arising  in  acoustics  and 
electromagnetics.  We  use  discrete  Fourier  techniques  (see  Richtmyer  and  Morton 
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[20])  to  analyze  the  difference  equations.  In  the  case  of  the  eigenvalue  problem,  we 
work  with  a  finite  domain  and  homogeneous  Dirichlet  boundary  conditions,  and 
for  the  Helmholtz  equation  we  perform  dispersion  analysis  on  infinite  domains, 
and  consider  a  boundary-value  problem  on  a  finite  domain. 

The  basis  of  comparison  used  throughout  this  paper  is  the  number  of  degrees  of 
freedom  in  the  discrete  model,  which  turns  out  to  be  equivalent  to  the  bandwidth 
of  the  corresponding  matrix  problem.  There  is  some  precedent  for  this  basis  of 
comparison,  namely,  it  was  used  by  Kwok,  Moser  and  Jimenez  [17]  in  studies  of 
B-splines,  finite  element,  and  collocation  methods  for  advective  and  diffusive  pro¬ 
cesses  and  they  presented  their  rationale  for  selecting  it.  Nevertheless,  one  can  still 
take  issue  with  it,  primarily,  in  our  opinion,  because  it  leads  to  significantly  differ¬ 
ent  numbers  of  quadrature  points  for  smooth  and  C°  basis  functions.  However,  it 
may  be  said  that,  for  the  smooth  case,  optimal  rules  are  not  yet  known  and  once 
they  are,  a  more  valid  comparison  of  cost  will  be  able  to  be  made.  In  the  mean¬ 
time,  we  will  use  the  number  of  degrees-of-freedom  as  a  basis  for  comparison,  but 
recognize  that  the  issue  is  more  complex. 

In  Section  2  we  briefly  review  the  problems  under  consideration.  In  Section  3  we 
recall  the  basis  of  isogeometric  analysis,  B-splines  and  NURBS.  We  describe  the 
different  geometric  constructions  which  lead  to  linear  and  nonlinear  parameteriza- 
tions  of  the  problem  domain.  This  has  important  consequences  in  vibration  analysis 
(see  [7]).  In  Section  4  we  begin  our  investigation  in  the  context  of  one-dimensional 
problems.  We  calculate  the  discrete  spectrum  of  the  eigenproblem,  and  the  dis¬ 
persion  properties  of  the  discrete  approximation  to  the  Helmholtz  problem  using 
complex  wave-number  analysis  [22].  After  presenting  the  details  for  linear  ele¬ 
ments,  we  state  a  “duality  principle”,  which  enables  us  to  map  results  of  spectral 
analysis  to  dispersion  analysis,  and  vice  versa.  Throughout  the  paper  we  invoke 
the  duality  principle  to  simplify  derivations.  Nevertheless,  there  are  subtle  differ¬ 
ences  between  spectrum  and  dispersion  analysis  that  need  to  be  noted,  namely, 
the  possible  existence  of  “outlier  frequencies”  [7]  in  spectrum  analysis,  and  the 
existence  of  complex  wave-numbers  leading  to  spurious  evanescent  waves  in  dis¬ 
persion  analysis.  These  phenomena  occur  for  higher-order  discretizations  and  so 
we  investigate  quadratics  in  some  detail,  and  sketch  what  happens  in  cubic  and 
higher-order  cases.  We  calculate  the  “stopping  bands”  for  classical  finite  elements, 
first  identified  by  Thompson  and  Pinsky  [22],  and  show  that  B- spline  s/NURBS  do 
not  engender  stopping  bands.  However,  they  produce  spurious  roots  corresponding 
to  evanescent  waves.  These  are  strongly  attenuated  and  do  not  seem  to  show  them¬ 
selves  in  numerical  calculations.  We  then  proceed  in  Section  5  to  a  two-dimensional 
model  problem  that  we  analyze  with  bilinear  elements.  This  problem  gives  us  the 
opportunity  to  explain  the  oscillations  in  frequency  errors  produced  in  numerical 
studies.  Based  on  results  in  Sections  3  and  4  we  are  able  to  confidently  use  numer¬ 
ics  to  calculate  invariant  analytical  spectra  for  classical  p-method  finite  elements 
and  NURBS.  This  comparison  is  quite  startling.  The  higher-order  p-elements  give 
rise  to  so-called  “optical  branches”  to  spectra,  which  have  no  approximation  prop- 
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erties.  It  is  well  known  that  the  upper  part  of  discrete  frequency  spectra  are  very 
inaccurate,  but  what  seems  to  be  a  completely  new  observation  is  the  errors  diverge 
with  p.  On  the  other  hand  there  are  no  optical  modes  with  NURBS  (at  least  when 
an  appropriate  “nonlinear”  parameterization  of  the  geometry  is  used  [7])  and  the 
spectral  errors  converge  with  p.  The  results  are  strikingly  different  (see  Figure  28) 
and  seem  to  register  a  significant  advantage  for  NURBS.  We  conjecture  that  these 
results  may  at  least  partially  explain  why  classical  higher-order  finite  elements  have 
not  been  widely  adopted  in  problems  for  which  the  upper  part  of  the  discrete  spec¬ 
trum  participates  in  a  significant  way,  such  as,  for  example,  impact  problems  and 
turbulence.  We  also  examine  eigenvectors  in  one  dimension  and  two-dimensional 
spectra  for  higher-order  approximations.  This  is  followed  by  studies  of  frequency 
response  spectra  and  wave  propagation  in  a  one-dimensional  rod.  In  all  cases,  we 
find  NURBS  outperform  standard  finite  elements.  Our  last  study  is  an  initiatory  one 
into  the  effects  of  reduced  numerical  quadrature.  We  find  that  reducing  quadrature, 
below  exact,  for  p-method  finite  elements,  has  deleterious  consequences.  Now,  not 
only  spectra  diverge  with  p,  but  for  a  fixed  p,  with  only  one  point  less  than  for  the 
exact  rule,  they  diverge  with  mesh  refinement  (i.e.,  //-refinement).  The  situation  is 
much  worse  than  for  the  exactly  integrated  case,  which  is  certainly  not  good.  It  is 
also  an  indication  that  the  Gauss  rules  are  indeed  optimal,  as  fewer  points  are  disas¬ 
trous.  On  the  other  hand,  reducing  quadrature  for  NURBS  does  not  have  significant 
effect.  Roughly  speaking,  the  Gauss  rule  with  approximately  half  as  many  Gauss 
points  as  required  for  exact  integration  provides  very  acceptable  results  with  no 
significant  degradation.  This  indicates  that  the  Gauss  rules  are  far  from  optimal  for 
NURBS.  In  fact,  it  is  conceptually  clear  that  the  Gauss  rules  are  not  the  answer  be¬ 
cause  they  do  not  acknowledge  in  any  way  the  precise  order  of  continuity  at  knots. 
So,  at  present,  application  of  the  Gauss  rules  between  knots  is  simple  and  effec¬ 
tive,  but  clearly  very  inefficient.  We  anticipate  that  optimal  rules  will  eventually  be 
developed  for  NURBS  and,  at  that  time,  we  will  be  able  to  make  more  equitable 
comparisons  of  cost.  We  draw  conclusions  in  Section  7. 


2  Structural  vibrations  and  wave  propagation 


In  this  section  we  briefly  recall  the  main  equations  of  structural  vibrations  and  of 
wave  propagation;  for  elaboration,  see  Chopra  [4],  Clough  and  Penzien  [5],  and 
Hughes  [14]  for  structural  vibrations;  Thompson  and  Pinsky  [23],  and  Thompson 
[24]  for  wave  propagation  (note  that  in  [23,  24]  particular  emphasis  is  on  acoustics). 
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2.1  Structural  vibrations:  natural  frequencies  and  modes 


Given  a  linear  (oo-dimensional)  structural  system,  the  undamped,  unforced  equa¬ 
tions  of  motion,  which  govern  free  vibrations,  are 

d2u 

M^  +  1C  u  =  0,  (1) 

where  M  and  /C  are,  respectively,  the  mass  and  stiffness  operators,  and  u  =  u(t ,  x) 
is  the  displacement. 

The  nth  normal  mode  <pn  and  its  frequency  t on  are  obtained  from  the  following 
eigenvalue  problem 

=  U2nM<(>n. 

We  remark  that  the  normal  modes  form  a  basis  in  space. 


Then,  we  can  separate  variables  as 

u(t,x)  = 

n 


and,  using  equation  (1),  obtain 


d  2  Unit) 
d  t2 


+  u2nun{t)  =  0. 


Each  mode  coefficient  un  oscillates  at  a  frequency  ujn,  and  we  can  write 


un  =  C_e~'lUnt  +  C+e'Unt 


After  discretization,  the  following  discrete  equations  of  motion  are  obtained 

d2uh 

M—  +  Kuh  =  0,  (2) 

where  M  and  K  are,  respectively,  the  finite-dimensional  consistent  mass  and  stiff¬ 
ness  matrices,  and  uh  =  uh(t.:  x)  is  the  discrete  displacement  vector. 

Analogously  to  the  continuum  case,  the  discrete  normal  modes  f>hn  and  the  frequen¬ 
cies  are  obtained  from  the  eigenproblem 

K  =  (w;)2M  <*>;,  (3) 

and  separating  variables  we  get 

uh(t,x)  =  E«nW0!l(x). 

n 
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with  oscillating  at  a  frequency  that  is, 


=  C_e 


-iu„t 


+  C4 


Aoji H 


The  nth  discrete  normal  mode  4>hn  is  an  approximation  of  the  nth  exact  normal  mode 
4>n,  for  n  —  1, . . . ,  TV,  being  N  the  total  number  of  degrees -of-freedom. 

The  corresponding  discrete  and  exact  frequencies  are  of  course  different  (see,  e.g., 
Figure  1). 


0  2  4  6  8  10  12  14  16  18  20 

mode  number 


Fig.  1.  Exact  and  discrete  natural  frequencies  for  the  one-dimensional  model  problem  of 
free  vibration  of  an  elastic  rod  with  homogeneous  Dirichlet  boundary  conditions.  The  dis¬ 
crete  method  is  based  on  linear  finite  elements. 


A  fundamental  question  is,  how  close  are  the  discrete  frequencies  to  the  continuous 
ones?  In  other  words,  how  well  does  the  discrete  spectrum  approximate  the  exact 
spectrum? 


2.2  Wave  propagation:  the  Helmholtz  equation 


The  classical  equation  governing  wave  propagation  is 


V2m 


1  d2w 
c2  dt2 


0, 


(4) 


where  c  is  the  wave  propagation  speed.  Particular  solutions  of  (4)  are  plane  waves 
of  frequency  cu  traveling  in  the  direction  n  at  a  speed  c,  which  can  be  expressed  as 
the  time-harmonic  wave  train 


w(x,t)  =  Re(Ael{-k™-ut)), 


(5) 
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where  k  =  uj/c  is  the  wave-number,  u  is  the  angular  frequency,  and  /I  is  a  complex 
number.  The  wavelength  (with  units  of  length)  is  defined  by  A  =  2n /k,  while  the 
dual  measure  of  period  (with  units  of  time)  is  defined  by  T  =  2tv/cu. 

Assuming  time-harmonic  solutions,  that  is,  with  abuse  of  notation,  u(t,  x)  =  elu>tu(x.), 
the  linear  wave  equation  (4)  reduces  to  the  Helmholtz  equation 

V2w  +  k2u  =  0,  (6) 

whose  solutions  in  Mn  are  linear  combinations  of  plane  waves  in  space  w(x)  = 
eifcn  x  After  discretization,  equation  (6)  gives  rise  to 

(K  -  k2M)uh  =  0.  (7) 

The  numerical  solution  of  the  above  equation  is  a  linear  combination  of  plane  waves 
having  numerical  wave-number  kh,  where,  in  general,  kh  ^  k. 

Thus,  discrete  and  exact  waves  have  different  wavelengths,  27t jkh  and  2i i/k  (see 
Figure  2). 


Fig.  2.  Different  exact  and  numerical  wave-numbers  produce  waves  with  different  wave¬ 
lengths. 

The  fundamental  issue,  which  is  addressed  by  dispersion  analysis,  is  to  determine 
the  dispersion  of  a  numerical  method,  that  is,  how  close  the  discrete  wave-number 
kh  is  to  its  continuous  counterpart  k. 


3  NURBS-based  isogeometric  analysis 


Non-Uniform  Rational  B-splines  (NURBS)  are  a  standard  tool  for  describing  and 
modeling  curves  and  surfaces  in  computer  aided  design  and  computer  graphics 
(see  Piegl  and  Tiller  [18]  and  Rogers  [21]  for  an  extensive  description  of  these 
functions  and  their  properties).  In  this  work,  we  use  NURBS  as  an  analysis  tool, 
which  is  referred  to  as  “isogeometric  analysis”  by  Hughes  et  al.  [15].  The  aim  of 
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this  section  is  to  present  a  brief  overview  of  features  and  properties  of  NURBS- 
based  isogeometric  analysis  for  ID  and  2D  problems.  We  will  utilize  NURBS  in 
our  study  of  the  problems  introduced  in  the  previous  sections.  The  section  starts 
with  a  short  description  of  B-splines  and  NURBS. 


3.1  B-splines  and  NURBS 


B-splines  in  the  plane  are  piecewise  polynomial  curves  composed  of  linear  combi¬ 
nations  of  B-spline  basis  functions.  The  coefficients  are  points  in  the  plane,  referred 
to  as  control  points. 


A  knot  vector  is  a  set  of  non-decreasing  real  numbers  representing  coordinates  in 
the  parametric  space  of  the  curve 

{£l  =  0)  •••!  in+p+l  =  1},  (8) 

where  p  is  the  order  of  the  B-spline  and  n  is  the  number  of  basis  functions  (and 
control  points)  necessary  to  describe  it.  The  interval  [£i,  £n+P+i]  is  called  a  patch.  A 
knot  vector  is  said  to  be  uniform  if  its  knots  are  uniformly-spaced  and  non-uniform 
otherwise.  Moreover,  a  knot  vector  is  said  to  be  open  if  its  first  and  last  knots 
are  repeated  p  +  1  times.  In  what  follows,  we  always  employ  open  knot  vectors. 
Basis  functions  formed  from  open  knot  vectors  are  interpolatory  at  the  ends  of  the 
parametric  interval  [0, 1]  but  are  not,  in  general,  interpolatory  at  interior  knots. 


Given  a  knot  vector,  univariate  B-spline  basis  functions  are  defined  recursively 
starting  with  p  =  0  (piecewise  constants) 


N,.  o(£) 


1  if  &  <  5  <  f(+i 

0  otherwise. 


(9) 


For  p  >  1  : 


N,A()  =  +  f+r+1_  f  AW,;;)-  do) 

S  i+p  Si  si+p+l  si+l 

In  Figure  3  we  present  an  example  consisting  of  n  =  9  cubic  basis  functions  gen¬ 
erated  from  the  open  knot  vector  (0,  0,  0,  0, 1/6, 1/3, 1/2,  2/3, 5/6, 1, 1, 1, 1}. 

If  internal  knots  are  not  repeated,  B-spline  basis  functions  are  -continuous.  If 
a  knot  has  multiplicity  k,  the  basis  is  C'p_fc-continuous  at  that  knot.  In  particular, 
when  a  knot  has  multiplicity  p,  the  basis  is  C°  and  interpolates  the  control  point  at 
that  location. 

By  means  of  tensor  products,  a  B-spline  region  can  be  constructed  starting  from 
knot  vectors  (£i  =  0, ...,  £n+P+i  =  1}  and  {rji  =  0,  ...,r]m+q+i  =  1},  and  annxm 
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Fig.  3.  Cubic  basis  functions  formed  from  the  open  knot  vector 

{0, 0, 0, 0, 1/6, 1/3, 1/2,  2/3,  5/6, 1, 1, 1, 1}. 

net  of  control  points  B?  J.  Two-dimensional  basis  functions  NitP  and  Mj  q  (with 
i  =  1, n  and  j  =  1, m)  of  order  p  and  q,  respectively,  are  defined  from  the 
knot  vectors,  and  the  B-spline  region  is  the  image  of  the  map  S  :  [0, 1]  x  [0, 1]  — >  Q 
given  by 

n  m 

SK.ij)  =  EE  "fcp  («)%,(»)By.  (ID 

i= 1  j=l 

The  two-dimensional  parametric  space  is  the  domain  [0, 1]  x  [0, 1].  Observe  that  the 
two  knot  vectors  (£i  =  0,  ...,£n+p+i  =  1}  and  {771  =  0,  ...,rjm+q+1  =  1}  generate 
in  a  natural  way  a  mesh  of  rectangular  elements  in  the  parametric  space. 

A  rational  B-spline  in  R2  is  the  projection  onto  two-dimensional  physical  space  of  a 
polynomial  B-spline  defined  in  three-dimensional  homogeneous  coordinate  space. 
For  a  complete  discussion  of  these  space  projections,  see  [6,  10]  and  references 
therein.  In  this  way,  a  great  variety  of  geometrical  entities  can  be  constructed  and, 
in  particular,  all  conic  sections  in  physical  space  can  be  obtained  exactly.  The  pro¬ 
jective  transformation  of  a  B-spline  curve  yields  a  rational  polynomial  curve.  Note 
that  when  we  refer  to  the  “order”  of  a  NURBS  curve,  we  mean  the  order  of  the 
polynomial  curve  from  which  the  rational  curve  was  generated. 

To  obtain  a  NURBS  curve  in  R2,  we  start  from  a  set  6  l3  (i  =  1, ... ,  n )  of 
control  points  (“projective  points”)  for  a  B-spline  curve  in  R3  with  knot  vector  H. 
Then  the  control  points  for  the  NURBS  curve  are 

/"DltA 

(B;)j  =  - — — ,  j  =  1,  2  (12) 

Wi 


where  (Bj)j  is  the  jth  component  of  the  vector  B;  and  u;,;  =  (B'';)3  is  referred  to  as 
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the  ith  weight.  The  NURBS  basis  functions  of  order  p  are  then  defined  as 


m) 


Ni)P(£)wi 

TJUNiJ^)wi 


(13) 


The  NURBS  curve  is  defined  by 


C(0  =  E^(0B,  (14) 

i= 1 


Analogously  to  B-splines,  NURBS  basis  functions  on  the  two-dimensional  para¬ 
metric  space  [0, 1]  x  [0, 1]  are  defined  as 


'/)  = 


Ni,p(€)Mjjq(r))wij 


that  the 

functions  are  the  same  as  for  B-splines. 


(15) 


where  Wij  =  (B™-)3.  Observe  that  the  continuity  and  support  of  NURBS  basis 


NURBS  regions,  similarly  to  B-spline  regions,  are  defined  in  terms  of  the  basis 
functions  (15).  In  particular  we  assume  from  now  on  that  the  physical  domain  Q 
is  a  NURBS  region  associated  with  the  n  x  m  net  of  control  points  By,  and  we 
introduce  the  geometrical  map  F  :  [0, 1]  x  [0, 1]  — ^  0  given  by 

n  m 

*= 1  3= 1 


3.2  Isogeometric  Analysis 


The  image  of  the  elements  in  the  parametric  space  are  elements  in  the  physical 
space.  The  physical  mesh  is  therefore 

%  =  {F((&,  &+i)  x  (rjj,  rjj+i)),  with  i  =  1, . . . ,  n  +  p,  j  =  1, . . . ,  m  +  q}  . 

(17) 

We  denote  by  h  the  mesh-size,  that  is,  the  maximum  diameter  of  the  elements  of 

rh. 

Following  the  isoparametric  approach,  the  space  of  NURBS  functions  on  O  is  de¬ 
fined  as  the  span  of  the  push-forward  of  the  basis  functions  (15) 

Vh  =  span{i?fj  o  (18) 
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3.2.1  Main  features 


In  the  following  we  present  a  summary  of  the  main  features  of  isogeometric  analy¬ 
sis.  The  interested  reader  may  find  more  details  and  applications  in  [2,  7,  8,  15,  19]. 

•  A  mesh  for  a  NURBS  patch  is  defined  by  the  product  of  knot  vectors. 

•  Knot  spans  subdivide  the  domain  into  “elements”. 

•  The  support  of  each  basis  function  consists  of  a  small  number  of  elements. 

•  The  control  points  associated  with  the  basis  functions  define  the  geometry. 

•  The  isoparametric  concept  is  invoked,  that  is,  the  unknown  variables  are  repre¬ 
sented  in  terms  of  the  basis  functions  which  define  the  geometry.  The  coefficients 
of  the  basis  functions  are  the  degrees-of-freedom,  or  control  variables. 

•  Three  different  mesh  refinement  strategies  are  possible:  analogues  of  classical  h- 
refinement  (by  knot  insertion)  and  /^-refinement  (by  order  elevation  of  the  basis 
functions),  and  a  new  possibility  referred  to  as  /.'-refinement,  which  increases 
smoothness  in  addition  to  order. 

•  The  element  arrays  constructed  from  isoparametric  NURBS  can  be  assembled 
into  global  arrays  in  the  same  way  as  finite  elements  (see  Hughes  [14],  chapter 
2). 

•  Dirichlet  boundary  conditions  are  applied  to  the  control  variables,  in  the  same 
way  as  in  finite  elements.  Neumann  boundary  conditions  are  satisfied  naturally 
as  in  standard  finite  element  formulations  (see  Hughes  [14],  chapters  1  and  2). 

Finally,  it  is  important  to  remark  that  in  structural  analysis  NURBS  elements  rep¬ 
resent  all  rigid  body  motions  and  constant  strain  states  exactly  (see  Hughes  [14]). 
Consequently,  structures  assembled  from  compatible  NURBS  elements  pass  stan¬ 
dard  “patch  tests”  (see  Hughes  [14],  chapters  3  and  4,  for  a  description  of  patch 
tests). 


3.3  Linear  and  nonlinear  parameterizations 


When  dealing  with  NURBS,  an  important  issue  is  the  choice  of  the  parameteriza¬ 
tion  to  be  used.  Take  as  an  example  a  ID  domain:  the  simplest  (and  more  natural) 
option  is  to  employ  a  linear  parameterization,  but  in  some  situations  a  nonlinear 
choice  can  be  more  suitable. 

The  isogeometric  procedure  originally  proposed  by  Hughes  et  al.  [15]  is  based  on 
a  distribution  of  control  points  which  leads  to  a  linear  parameterization  (i.e.,  con¬ 
stant  Jacobian  determinant),  but  in  Cottrell  et  al.  [7]  it  has  been  shown  that  when 
studying  structural  vibrations  a  nonlinear  parameterization,  such  that  the  control 
points  are  uniformly  spaced,  gives  better  results.  In  Figure  4,  we  show  the  ID  dis¬ 
tribution  of  21  control  points  obtained  for  the  two  cases  using  cubic  NURBS  (top), 
along  with  plots  of  the  corresponding  parameterization  x  =  x(f  )  and  Jacobian 
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d  3c(£^) 

J(£)  =  g  £  (bottom).  Subsequently,  we  will  refer  to  this  choice,  in  which  con¬ 
trol  points  are  uniformly  distributed,  as  “nonlinear  parameterization”,  in  contrast 
with  the  linear  one. 


0  .1  .2  .3  .4  .5  .6  .7  .8  .9  1 

X 


Fig.  4.  ID  case:  linear  versus  nonlinear  parameterization  determined  by  uniformly-spaced 
control  points  (cubic  NURBS,  21  control  points).  Top:  distribution  of  control  points;  dots 
correspond  to  linear  parameterization  control  points  and  asterisks  to  uniformly-spaced  con¬ 
trol  points.  Bottom:  plot  of  the  parameterization  (left)  and  of  its  Jacobian  (right)  for  the  two 
cases. 


Finally,  referring  to  the  2D  case,  we  present  in  Figure  5  the  example  of  a  control 
net  and  mesh  (i.e.,  the  physical  representation  of  the  elements)  of  a  square  physical 
domain,  obtained  using  p  =  q  =  A  and  11x11  control  points  for  both  the  linear 
and  the  nonlinear  parameterizations. 


3.4  k-method  and  p-method 


We  conclude  this  section  on  isogeometric  analysis  by  briefly  pointing  out  what 
we  mean  in  this  paper  by  the  terms  “k-method”  and  “p-method”.  Referring  to  the 
already  cited  k-  and  p-refinement  strategies,  we  define  the  k-method  as  the  analysis 
method  exploiting  the  full  continuity  across  the  elements  allowed  by  NURBS  basis 
functions  (i.e.,  Cp~l  for  a  degree  p  NURBS).  In  the  following  we  will  simply  label 
this  method  as  “NURBS”.  Instead,  we  define  the  p-method  as  the  analysis  method 
where  only  (7° -continuity  is  enforced  across  elements  (this  can  be  obtained  with 
isogeometric  analysis  by  repeating  the  knots  of  a  degree  p  NURBS  p  —  1  times). 
This  approach,  used  in  combination  with  a  linear  parameterization,  is  equivalent  to 
classical  hp- finite  element  methods,  and  in  the  following  we  will  simply  label  it  as 
“FEM”. 
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Fig.  5.  2D  case:  linear  versus  nonlinear  parameterization  determined  by  uniformly-spaced 
control  points  (p  =  q  =  4,  11  x  11  control  points).  Top:  control  net  (left)  and  mesh  (right) 
obtained  employing  the  linear  parameterization,  both  plotted  on  the  physical  domain.  Mid¬ 
dle:  mesh  on  the  parent  domain.  Bottom:  control  net  (left)  and  mesh  (right)  obtained  em¬ 
ploying  the  nonlinear  parameterization,  both  plotted  on  the  physical  domain. 

4  Analytical  study  in  one  dimension 


In  this  section,  we  carry  out  some  analytical  computations  for  finding  the  discrete 
spectrum  for  structural  vibrations  (spectrum  analysis)  and  the  dispersion  relation 
for  wave  propagation  (dispersion  analysis),  and  we  discuss  the  similarity  between 
the  two  frameworks.  We  first  deal  with  the  case  of  an  approximation  with  linear 
elements,  for  which  k-  and  p- methods  coincide.  Then,  we  discuss  the  extension  of 
the  results  to  higher  order  approximations,  for  both  methods.  This  section  is  partly 
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based  on  [1,  12-14,  16,  22]. 


4.1  Linear  approximation 

4.1.1  Spectrum  analysis 

In  this  section,  we  solve  the  generalized  discrete  eigenproblem  (3)  associated  to 
a  linear  approximation  on  the  one-dimensional  domain  (0,  L ),  and  study  the  error 
between  discrete  and  exact  solutions.  We  employ  a  uniform  mesh  0  =  xq  <  x\  < 
. . .  <  xa  <  ■  ■  ■  <  xN+i  =  L,  where  the  number  of  elements  is  nei  =  N  +  1  and 
the  mesh-size  is  h  =  L/nei. 

Considering  homogeneous  Dirichlet  (fixed-fixed)  boundary  conditions,  the  eigen¬ 
problem  (3)  can  be  written  as 

1  h(coh)2 

yi^A-i  ~  Z&A  +  rf’A+i)  H - - — ((^a-i  +  ^a+^a+i)  —  0,  A  —  1, ...  ,N,  (19) 

h,  o 

0o  =  0iv+i  =  0,  (20) 

where  N  is  the  total  number  of  degrees-of-freedom,  and  0a  =  4>h{ xa )  is  the  nodal 
value  of  the  discrete  normal  mode  at  node  xa-  Equation  (19)  constitutes  a  lin¬ 
ear  homogeneous  recurrence  relation  of  order  2,  whose  solutions  (ignoring,  for 
now,  the  boundary  conditions  (20))  are  linear  combinations  of  exponential  func¬ 
tions  4>a  =  ( p\)A  and  (j)A  =  (p2)A,  where  p\  and  p2  are  the  distinct  roots  of  the 
characteristic  polynomial 


(1  -  2p  +  p2)  +  ^-^-(1  +  4 p  +  p2)  =  0.  (21) 

o 

Actually,  (21)  admits  distinct  roots  when  a )hh  ^  0,  a/T2;  for  ujhh  =  0,  (21)  admits 
the  double  root  p  —  1  (in  this  case,  solutions  of  (19)  are  combinations  of  0^  =  1 
and  (f>A  =  A,  that  is,  the  affine  functions),  while  for  ujhh  =  \/l2  there  is  a  double 
root  p  —  —  1  (and  solutions  of  (19)  are  combinations  of  0^  =  (— 1)A  and  0^  = 
A{— 1)A).  Observe  that,  in  general,  p2  =  pi1. 

For  the  purpose  of  spectrum  analysis,  we  are  interested  in  0  <  u>hh  <  Vl2,  which 
we  assume  for  the  remainder  of  this  section.  In  this  case,  p12  are  complex  conjugate 
(we  assume  Im(pi)  >  0)  and  of  unit  modulus.  Moreover,  in  order  to  compare  the 
discrete  spectrum  to  the  exact  spectrum,  it  is  useful  to  represent  the  solutions  of 
(19)  as  linear  combinations  of  e±lAujh  (that  is,  4>a  =  C-e~lAujh  +  C+elAujh),  by 
introducing  uj  such  that  eluh  =  p\.  With  this  hypothesis,  to  is  real  and,  because  of 
periodicity,  we  restrict  to  0  <  c oh  <  ir.  Using  this  representation  in  (21)  and  using 
the  identity  2  cos(o)  =  eia  +  e~ia,  after  simple  computations  the  relation  between 
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uh  and  uhhis  obtained: 


(uhh)2 

6 


(2  +  cos  (uh))  —  (1 


cos  (uh))  =  0. 


Solving  for  uhh  >  0,  we  get 


(22) 


c uhh 


1  —  cos  (uh) 

2  +  cos  (uh) 


(23) 


Taking  now  into  account  the  boundary  conditions,  non-null  solutions  4>a  of  (19)- 
(20)  exist  when  a;  —  n/L,2n/L, . . . ,  Nn/L.  Indeed,  for  uh  =  nnh/L  =  nn/(N+ 
1),  and  C_  =  —  C+,  the  solution 


4>a 


e+iAnir/(N+l ) 

c - 


_  g— iAnn/(N+l) 

2i 


C  sin 


Ann  \ 

ivTTy 


vanishes  when  A  =  QorA  =  NJ\-l. 


(24) 


Precisely,  (24)  is  the  nth  discrete  normal  mode,  associated  to  the  corresponding  nth 
discrete  natural  frequency  uh,  given  by  (23): 


N  +  l  1  —  cos(nn/(N  +  1)) 
L  \  2  +  cos(mr  /  (N  +  1)) 


(25) 


Observe  that  (25)  returns  the  frequencies  u>h  in  increasing  order  with  respect  to 
n.  Figure  6  shows  the  dimensionless  discrete  natural  frequencies  a Ah,  for  N  = 
9  degrees-of-freedom.  They  are  represented  by  points  lying  on  the  graph  of  c ohh 
versus  uh/n,  given  by  (23).  The  abscissa  uh/n  is  equivalent  to  the  scaled  wave- 
number  n/(N  +  1). 


As  is  known,  the  nth  discrete  mode  (j)A 
terpolant  of  the  nth  exact  mode  cf)(x)  = 

ojh  a ; 

is  u  =  nn / L.  The  quantity - l  —  — 

UJ 

natural  frequency.  The  plot  of 


=  Csm(Ann/(N  +  1))  is  the  nodal  in- 
C  sin(nnx/L),  whose  natural  frequency 

h  —  ui 

- represents  the  relative  error  for  the 

UJ 


UJ 


h 


UJ 


l 

uh \ 
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1  —  cos(uh ) 

2  +  cos(uhy 


(26) 


straightforwardly  derived  from  (23),  is  shown  in  Figure  10. 


4.1.2  Dispersion  analysis 

We  obtain  here  the  discrete  dispersion  relation  for  linear  approximation.  We  con¬ 
sider  the  Helmholtz  equation  (6)  on  the  infinite  domain  (line),  and  its  discretization 
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Fig.  6.  Analytically-computed  (discrete)  natural  frequencies  for  linear  approximation, 
N  =  9.  The  dimensionless  frequencies  lie  on  the  graph  (dotted  line)  of  relation  (23). 

(7)  on  the  numerical  grid  xA  =  hA,  AeZ.  The  resulting  stencil  equation  is 


As  we  described  in  Section  2.2,  the  standard  dispersion  analysis  consists  of  com¬ 
paring  the  wave-numbers  of  the  exact  and  discrete  solutions.  We  recall  that  the 
exact  solutions  of  the  Helmholtz  equation  are  linear  combinations  of  u(x)  =  e±lkx. 
Also,  the  discrete  solutions,  that  is,  solutions  of  the  stencil  equation  (27),  are  com¬ 
binations  of  exponentials  (as  we  have  seen  in  the  previous  section  for  (19),  which 
is  analogous  to  (27)):  following  the  notation  which  is  common  in  the  context  of 
dispersion  analysis,  the  discrete  solutions  are  written  as 


uA  =  uh{xA)  =  C_e~ikhhA  +  C+eikhhA 


(28) 


denoting  by  kh  the  discrete  wave-number.  In  general,  kh  6  C,  and  is  uniquely 
determined  under  the  condition  0  <  Re(khh)  <  n.  Inserting  (28)  into  (27),  the 
relation  between  k  and  kh  is  obtained  (analogously  to  (22))  as 


(29) 


In  this  context,  one  is  usually  interested  in  solving  (29)  with  respect  to  khh.  The 
first  step  is 


(30) 
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where  it  is  easily  seen  that  (since  kh  >  0) 


6-2  {kh)2 
6  +  {kh)2 


<=>•  kh  <  \/l2. 


Then,  for  kh  <  a/12  the  discrete  wave-number  is  real  and  given  by 


khh  =  arccos 


/6  — 2  {kh)2\ 
\6+  {kh)2  J 


(31) 


when  kh  >  a/T2,  khh  can  be  obtained  as  in  (31),  but  arccos(-)  has  to  be  understood 
as  the  complex  arc-cosine  (see  [22]).  The  non-zero  imaginary  part  of  khh  produces 
an  amplitude  modulation  of  the  discrete  solutions  which  is,  clearly,  an  unphysical 
feature  of  the  numerical  solution.  However,  for  linear  elements,  it  happens  when 
reaching  the  resolution  limit ,  which  corresponds  to  the  largest  wave-number  that 
the  numerical  mesh  can  represent  (before  aliasing  occurs).  For  linear  elements,  the 
resolution  limit  is  khh  =  it. 

The  plots  of  Re{khh)  and  Im{khh)  versus  kh  are  shown  in  Figure  7  (observe  that 
kh  is  represented  on  the  ordinate).  The  amplitude  spectrum,  that  is,  \u/^f\uA+\  = 
eim(k  h)  versus  (assuming  ua  =  elk  hA),  is  presented  in  Figure  8. 


Fig.  7.  Analytically-computed  (discrete)  wave-number  for  linear  approximation. 

The  dispersion  error  {kh  —  k)/k  =  kh / k  —  1  is  typically  displayed  in  the  litera¬ 
ture  (e.g.,  see  [12])  by  plotting  the  quantity  kh  jk  versus  khh,  for  khh  real.  This  is 
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Fig.  8.  Amplitude  spectrum  for  linear  approximation, 
obtained  from  (29)  as 


kh 

T 


khh 


1  —  cos  (khh) 
\  2  +  cos  (khh) 


(32) 


and  is  shown  in  Figure  9. 


Fig.  9.  Discrete-to-exact  wave-number  ratio  for  linear  approximation. 
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Fig.  10.  Unified  dispersion  and  spectrum  analysis  for  linear  approximation. 

4.1.3  Duality  Principle 

Note  that  (32)  is  the  reciprocal  of  expression  (26),  written  in  terms  of  different 
quantities:  to  corresponds  to  kh  while  u>h  corresponds  to  k. 

Indeed,  it  is  clear  from  Section  4. 1.1-4. 1.2  that  spectrum  analysis  is  equivalent  to 
dispersion  analysis  in  the  regime  where  kh  is  real:  switching  from  one  field  to  the 
other  is  just  a  matter  of  exchanging  notation,  from  the  mathematical  viewpoint. 

From  now  on,  we  will  represent  the  dispersion  error  by  plotting  the  ratio  k/kh  ver¬ 
sus  khh.  While  this  is  not  common  in  literature,  it  is  suitable  for  unifying  dispersion 
and  spectrum  analysis  (see  Figure  10). 

Remark  1  Figure  10  can  be  obtained  numerically,  instead  of  analytically.  After 
numerically  computing  the  spectrum,  with  eigenvalues  sorted  in  increasing  order, 
and  then  the  frequencies  of  the  discrete  system,  each  discrete  frequency  is  divided  by 
the  corresponding  exact  frequency  rm  /  L.  This  gives  uh/co,  with  correct  association 
of  discrete  to  exact  modes.  In  this  case,  the  scaled  mode  number  n/(N  +  1 )  has  to 
be  represented  as  the  abscissa. 


4. 2  Higher  order  p-method 


We  have  seen  in  the  previous  section  that,  for  linear  elements,  spectrum  and  disper¬ 
sion  analysis  are  equivalent.  This  holds  for  the  p-method  with  higher  order  elements 
as  well,  though  the  analysis  becomes  more  technical.  We  discuss  first,  and  in  more 
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detail,  the  simpler  case  of  quadratic  elements. 

We  start  by  focusing  on  dispersion  analysis.  Therefore,  we  take  into  consideration 
the  Helmholtz  equation  on  the  infinite  line  and  its  discretization  by  quadratic  finite 
elements  on  the  uniform  grid  . . .  <  xa  <  %a+ 1/2  <  %a+ 1  <  ■  ■ where  xa  =  hA 
(for  A  G  Z)  are  the  element-endpoint  nodes  and  xa+i/2  =  h(A  +  1/2)  (for  A  G  Z) 
are  the  mid-point  nodes.  On  this  mesh,  we  consider  the  usual  nodal  basis,  depicted 
in  Figure  11.  The  corresponding  stencil  equation  is  different  for  element-endpoint 


Fig.  11.  Basis  functions  for  the  quadratic  p-method. 
degrees-of-freedom  and  bubble  (internal  to  element)  degrees-of-freedom:  one  has 

7rr{— uA-1  +  8ua~  1/2  —  14  rt  A  +  8Ua+ 1/2  —  Wjt+l) 

6h  (33) 

+  k2—(  —  UA-l  +  2UA-1/2  +  8ua  +  <2,Ua+1/2  ~  UA+ 1)  =  0,  V54  G  Z. 


and 

1 


h 


&ua  ~  IQua+i/2  +  8ma+i)  +  k"  —[2ua  +  16w A+1/2  +  211,4+1)  —  0,  VA  G  Z, 

(34) 

respectively.  One  could  look  for  a  solution  of  (33)— (34)  at  the  element-endpoint 
and  bubble  nodes  (see,  for  example,  [22,  §3.2.4]).  However,  a  simpler  and  more 
common  way  to  proceed  consists  of  calculating  the  bubble  degrees-of-freedom  as 


UA+l/2  ~ 


40  +  ( khf 


-(ua  +  Ua+i), 


8(10-  (kh)2j 

and  eliminating  them,  obtaining  a  system  of  equations  for  ua,  A  G  Z,  which  is 


(35) 


1 

3 h 


k2- 

30 


'30  +  2  (khf 
v  10  -  (kh)2 

5  (kh)' 


Ua-1  + 


-60  +  16  (khf 


ua 


'30  +  2  {khf 


v40-4  (khf 


U  A — 1  + 


10  -{khf  j 
/200  -  15  {khf 


20  -  2  {khy 


ua  + 


10-  (kh)2  J 
/  5  (kh)2 


ua+ 1 


v40-4  (khf 


ua+ 1 

(36) 


=  0, 
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Remark  2  The  bubble  elimination  is  not  possible  when  the  bubble  equation  (34) 
is  singular  for  ua+ 1/2-  VVfe  refer  to  this  situation  as  bubble  resonance.  It  occurs  for 

kh  =  /Lb.  (37) 


Observe  that  (36)  is  a  homogeneous  linear  recurrence  equation  of  order  2,  as  for 
the  linear  case  (27).  Then,  its  solutions  can  be  written  as 

uA  =  C-e~ikhhA  +  C+eikhhA,  \/A  e  Z.  (38) 


Substituting  (38)  into  (36),  one  obtains 


cos  (khh) 


3  khA  —  104  kh2  +  240 
khA  +  16  kh 2  +  240 


(39) 


Given  kh,  there  is  only  one  solution  khh  to  (39)  if  we  restrict  to  0  <  Re(khh)  <  7 r; 
this  is 


khh  =  arccos 


/  3  kh4  —  104  kh2  +  240  \ 

\  kh 4  +  16  kh2  +  240  )  ' 


(40) 


As  for  the  linear  case,  if  the  right  hand  side  of  (39)  is,  in  modulus,  smaller  than  or 
equal  to  1,  then  khh  is  real  (see  the  left  plot  of  Figure  12).  From  (40)  and  Figure 


Fig.  12.  Analytically-computed  (discrete)  wave-number  for  quadratic  p-method,  and  khh 
real.  Relation  (40)  is  plotted  on  the  left,  and  relation  (41)  on  the  right. 

12  (left),  it  is  seen  that  each  real  value  of  kfh  is  associated  with  two  values  of 
kh,  on  two  different  branches,  termed  acoustical  and  optical  (cf.  [3]).  This  means 
that  the  solution  at  the  endpoint-element  nodes  of  the  grid  is  the  same  for  the  two 
corresponding  kh’s;  however,  the  bubble  degrees-of-freedom  (given  by  (35))  are 
different  for  the  two  cases,  which  means  that  the  two  discrete  solutions  (at  element- 
endpoint  and  bubble  nodes)  are  different.  From  (40),  it  is  found  that  kh  G  [0,  vTO] 
for  the  acoustical  branch  and  kh  e  [Vl2,  a/60]  for  the  optical  branch. 
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Fig.  13.  Analytically-computed  (discrete)  wave-number  for  quadratic  p-method. 

A  better  representation  of  the  relation  khh  versus  kh  is  derived  from  (39)  as  fol¬ 
lows.  For  the  optical  branch,  that  is,  for  kh  >  Vl2,  thanks  to  the  even  parity  and 
periodicity  of  the  cosine  function,  we  represent  the  discrete  wave-number  in  the 
range  n  <  Re(khh )  <  27 r.  Then,  we  set 


khh 


arccos 


f  3  kh4  —  104  kh2  +  240  \ 

,  kh4  +  16  kh2  +  240  / 

3  kh4  —  104  kh2  +  240 


arccos 


UUA  i  1  R  UU 2 


for  kh  <  a/12, 
for  kh  >  \/l2. 


(41) 


This  results  in  a  one-to-one  monotone  relation  between  kh  and  the  real  values  of 
khh,  plotted  in  Figure  12  (right),  which  is  consistent  with  the  physical  expectation 
and  is  useful  in  view  of  spectrum  analysis.  Moreover,  khh  =  27t  plays  the  role 
of  resolution  limit  of  the  complete  numerical  grid  (considering  both  bubble  and 
element-endpoint  nodes). 

Allowing  complex  wave-numbers  in  (41),  Figure  13  is  obtained.  Notice  that  Im(khh) 
is  not  only  zero  for  kh  >  a/60  but  also  in  between  the  two  branches,  for  \/l0  < 
kh  <  a/I2.  This  interval  is  called  a  stopping  band  and  its  effect  on  the  numerical 
solution  will  be  discussed  in  Section  6.2.1.  The  amplitude  spectrum  is  shown  in 
Figure  14. 

We  now  turn  to  spectrum  analysis,  and  consider  the  eigenproblem  (3)  on  the  do¬ 
main  (0,  L ).  The  mesh  restricts  to  0  =  x0  <  . . .  <  xa-i/2  <  xa  <  %a+ 1/2  < 

. . .  <  xUel  =  L.  The  mesh- size  is  h  =  L/nei.  We  have  therefore  nej  +  1  element- 
endpoint  nodes,  including  x0  and  xrifd .  and  nc!  bubble  nodes.  Taking  into  account 
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Fig.  14.  Amplitude  spectrum  for  quadratic  p-method. 


the  homogeneous  Dirichlet  boundary  conditions,  there  are  N  =  2 nei  —  1  degrees- 
of-freedom.  Based  on  the  previous  study  of  the  Helmholtz  equation  (33)— (34),  we 
assume  ujhh  \/l0  and  perform  the  bubble  elimination,  leading  to  the  following 
equation  for  the  element-endpoint  degrees-of-freedom:  for  A  =  1, . . . ,  ne/  —  1 


(42) 


We  also  have  the  boundary  conditions  </>0  =  =  0.  Clearly,  (42)  is  the  coun¬ 

terpart  of  (36).  We  proceed  now  using  the  dispersion  analysis  results  of  the  present 
section,  invoking  the  duality  principle,  that  is,  the  correspondence  ujh  ^  k  and 
ijj  kh,  and  reasoning  as  for  spectrum  analysis  in  the  linear  case  (see  Section 
4.1.1). 

Normal  modes  at  element-endpoint  nodes  can  be  written  as  (j)A  =  C_e~VjjhA  + 
C+elujhA ;  the  boundary  condition  (f> o  =  0  determines  C_  =  —  C+,  while  (j>nel  =  0  if 
—  G  Z.  Observe  that  the  complex  values  of  ujh  are  not  of  interest  in  this  case.  The 
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relation  between  uhh  and  uh  is  (analogous  to  (39)) 


,  , ,  3  (c ohhf  -  104  (c ohh)2  +  240 

COS(uh)  =  — — — — - - — j — rr; - 

(uhh)4  +  16  ( uhh )2  +  240 


(43) 


The  natural  frequencies  are  obtained  solving  (43)  with  respect  to  uhh.  Unlike  the 
linear  case,  this  is  a  two-branch  relation:  there  are  two  roots  uhh  >  0  for  any  c oh  6 
M.  As  we  have  seen,  a  monotone  uhh  versus  uh  relation  is  obtained  representing 
the  two  branches  in  the  range  uh  E  [0,  ir\  and  uh  G  [ir,  27t]  ,  respectively.  Therefore, 
we  associate  to 

77.7T 

Ujh= — ,  n  =  1, . .  .nei  -  1,  (44) 

nei 

the  smallest  positive  root  of  (43),  obtaining  the  acoustical  branch,  and  we  associate 
to 

TL7T 

uh  =  — ,  n  =  nei  +  1, . . .  2nei  —  1  =  N]  (45) 

nei 

the  highest  root  of  (43),  obtaining  the  optical  branch.  These  roots  are  the  natural 
frequencies  that  can  be  obtained  by  bubble  elimination. 

The  frequency  uhh  =  vTO,  which  gives  bubble  resonance  (see  Remark  2)  has  to 
be  taken  into  consideration  as  well.  Indeed,  it  is  associated  with  the  normal  mode 

(t>A  =  0,  \/A  =  0,...,nei, 

<t>A+i/2  =  C{-  1)A  VA  =  0, . . .  ,nei  —  1. 

Since  uhh  =  \/l0  is  located  between  the  two  branches,  this  frequency  is  associated 
with  mode  number  n  =  nei-  Notice  that  with  this  choice,  all  the  normal  modes  at 
element  endpoints  are  given  by 

N+  -J  ,  A  =  0, 1, . .  .nei,  (47) 

n  being  the  mode  number.  Therefore,  (47)  interpolate  of  the  exact  modes  (at  ele¬ 
ment  endpoint  nodes). 

Eventually,  this  results  in  a  monotone  ordering  of  all  the  natural  frequencies  of  the 
discretized  system,  as  shown  in  Figure  15  (for  N  =  9).  The  abscissa  in  Figure  15 
is  uh/2TT,  which  corresponds  to  the  scaled  wave-number  n/  (N  +  1). 

The  numerical  error  in  the  calculation  of  natural  frequencies  is  visualized  by  the 
graph  of  uh/u  versus  uh.  As  for  the  linear  case,  it  is  the  same  graph  of  k/kh  versus 
khh,  which  reveals  the  dispersion  error.  The  unified  plot  is  shown  in  Figure  16. 

What  we  have  described  for  the  case  p  =  2  can  be  easily  generalized  to  higher 
order  cases  p  >  2,  in  one  space  dimension.  In  general  we  have  p  —  1  bubbles  per 
element,  and  then  p  branches  and  p  —  1  stopping  bands.  See,  for  example,  Figures 
17  and  18  which  refer  to  the  case  p  —  3.  In  particular,  assuming  the  resolution  limit 
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Fig.  15.  Analytically-computed  (discrete)  natural  frequencies  for  quadratic  p-method 
{N  =  9). 
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Fig.  16.  Unified  dispersion  and  spectrum  analysis  for  quadratic  p-method. 


to  be  khh  =  toh  =  pit,  one  can  derive  the  monotone  tohh  versus  a ih  relation  which 
is  useful  for  the  spectrum  representation.  There  are  also  p  —  1  natural  frequencies 
associated  to  bubble  resonance. 
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Fig.  17.  Analytically-computed  (discrete)  wave-number  for  cubic  p-method. 
4.3  Higher  order  k-method 


There  is  also  a  strict  relation  between  spectrum  and  dispersion  analysis  for  the 
/c-method  with  higher-order  elements.  Again,  for  the  sake  of  simplicity,  we  only 
discuss  in  detail  the  quadratic  approximation,  and  briefly  mention  the  extensions  to 
p  >  2  at  the  end  of  this  section. 

Considering  the  Helmholtz  equation  (7)  on  the  infinite  line,  we  denote  now  by 
xa  =  hA,  for  A  e  Z,  the  sequence  of  /i-spaced  control  points  (giving  a  linear 
parameterization  of  the  infinite  line);  the  stencil  equation  is  then 


~gr(uA- 2  +  2uA.l  —  Qua  +  2ua+i  +  uA+2) 

6h  h  (48) 

+  k 2  ^  (it a— 2  +  26-Uyi-i  +  66ua  +  2614,4+1  +  M4+2)  =  0,  VA  E  Z. 

Remember  that,  for  the  /r- method,  uA  denotes  the  coefficient  in  the  basis  expan¬ 
sion,  which  is  no  longer  interpolatory.  Another  major  difference  from  the  cases 
considered  in  the  previous  sections  is  that  (48)  is  a  homogeneous  recurrence  rela¬ 
tion  of  order  4.  Because  of  its  structure,  its  solutions  can  be  written  as  linear  com¬ 
binations  of  the  four  solutions  e±lkhflA  and  e±lkhhA,  where  kh  ^  kh  are  uniquely 
determined  under  the  assumption  0  <  Re(khh )  <  7T,  0  <  Re(khh)  <  7r  and 
Im(kh )  <  Im(kh )  <  0.  Notice  that  the  space  of  discrete  solutions  has  dimension 
4,  unlike  the  space  of  exact  solutions  (e±lkx)  which  has  dimension  2.  The  values 
of  khh  and  khh  can  be  obtained  from  kh  using  (48).  These  are  plotted  in  Figure 
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Fig.  18.  Amplitude  spectrum  for  cubic  p-method  (top)  and  detail  of  the  two  stopping  bands. 


19.  It  is  seen  that  kh  is  an  approximation  of  the  exact  wave-number  k  ( khh  ~  kh 
within  the  resolution  range  khh  e  [0, 7r]),  while  kh  is  a  numerical  wave-number, 
associated  with  spurious  evanescent  waves  of  the  form  ua  =  C(—l)AeZfIm('khh'>A. 

The  role  of  the  spurious  solutions  of  (48)  is  not  fully  clear.  These  solutions  are 
irrelevant  at  low  kh  (in  this  case,  the  a  priori  error  analysis  guarantees  the  accuracy 
of  the  numerical  solution)  while  they  could  affect  the  numerical  solution  at  high 
kh.  Nevertheless,  in  all  the  numerical  tests  we  have  performed  (see  Section  6.2), 
they  did  not  appear,  perhaps  because  they  are  so  strongly  attenuated.  See  Figure  20. 

Therefore,  for  the  purpose  of  the  dispersion  analysis,  we  only  consider  the  relation 
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Fig.  19.  Analytically-computed  (discrete)  wave-number  for  quadratic  /.'-method. 


Fig.  20.  Analytically-computed  spurious  wave-number  for  quadratic  /.  -method. 
k/kh  versus  khh.  After  simple  computations  this  is  obtained  from  (48)  as 


k 


20(2  —  cos  (khh)  —  cos  (khh)2) 


(49) 


kh  khh,\  16  +  13  cos  (khh)  +  cos  (khh)2  1 

and  plotted  in  Figure  21. 

Let  us  discuss  now  how  the  results  above  are  related  to  spectrum  analysis.  Consider 
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Fig.  21.  Unified  dispersion  and  spectrum  analysis,  for  quadratic  /c-method. 

the  eigenvalue  problem  (3)  on  the  interval  (0,  L ),  on  which  we  introduce  N+2  con¬ 
trol  points  associated  to  a  linear  parameterization.  The  control  points  are  uniformly 
spaced,  at  distance  h,  only  in  the  interior  portion  of  the  domain,  while  they  get 
closer  to  each  other  at  the  endpoints  of  the  domain  (0,  L ),  as  shown  in  Figure  4. 
On  the  other  hand,  the  space  of  discrete  functions  is  made  of  piecewise  quadratic 
polynomials  on  a  uniform  mesh  with  knot  spacing  h 1 ,  with  global  C 1  regularity 
(see  Figure  22).  In  this  cases,  the  discrete  problem  (3)  reads 


I 


(50) 


VA  =  3, . . . ,  N  -  2, 

7rr(0A-2  +  20,4-1  —  6  (ft  a  +  Z^a+i  +  (t)A+ 2) 

bn 


(51) 


1  Note,  h  here  corresponds  to  the  knot  spacing  in  physical  space.  Previously,  when  analyz¬ 
ing  quadratic  finite  elements,  we  adopted  the  usual  convention  that  h  represents  the  element 
length.  Consequently,  the  h’ s  in  these  two  cases  differ  by  a  factor  of  2. 
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7rr(0iV-3  +  20at_2  —  607V-1  +  4>n  +  20tv+i) 

bn 

H - 1  on  (0-ZV-3  +  260AT-2  +  660AT-1  +  25  07V  +  20at_|_i)  =  0, 

<  x  120  (52) 

7rr(4>N- 2  +  0JV-1  —  807V  +  607V+l) 

bn 

h(uh)2 

H -  {4>n-2  +  250AT-1  +  4007V  +  140tv+i)  =  0 

IzU 

along  with  the  Dirichlet  boundary  conditions 

00  =  0iV+l  =  0.  (53) 

Substituting  (53)  into  (50)-(52),  this  is  an  iV-dimensional  generalized  eigenvalue 
problem. 

The  equations  (51)  for  the  interior  degrees -of-freedom  correspond  to  the  stencil 
(48)  we  previously  considered  in  the  dispersion  analysis.  However,  the  boundary 
equations  (50)  and  (52)  are  different,  due  to  the  different  shape  of  the  boundary 
basis  functions  (see  Figure  22).  The  relation  between  the  solutions  of  (48)  and  the 


Fig.  22.  Basis  functions  for  the  quadratic  fc-method 

ones  of  (50)-(53),  which  is  not  trivial,  is  addressed  in  the  following  discussion. 
First,  we  deal  with  (51)  and  use  the  previous  study  of  the  Helmholtz  equation  (in¬ 
voking  the  duality  principle  and  change  of  notation  tuh  k,  to  kh  and  Co  <-►  kh), 
to  infer  that  the  solutions  of  (51)  are 

0A  =  C+eluhA  +  C.e~iujhA  +  C+j*hA  +  C-e~iCbhA,  VA  =  1,...,N.  (54) 


Inserting  this  expression  into  the  boundary  equations  (50)  and  (52),  and  imposing 
the  boundary  conditions  (53),  one  determines  the  four  coefficients  C_,  C+,  C_, 
and  C+  in  order  that  (54)  is  the  solution  of  the  entire  system  (50)-(53).  The  trivial 
solution  corresponds  to  C_  =  C+  =  C-  =  C+  =  0.  The  normal  modes  are  instead 
the  non-zero  solutions,  which  exist  for  suitable  uohh  (giving  the  discrete  natural 
frequencies)  and  corresponding  uoh  and  uh.  Precisely,  the  nth  normal  mode  turns 
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(55) 


out  to  be  2 


that  is,  uh 
given  by 


|  0a  =  C  sin  (^rmA  ,  VA  =  1, . . . ,  N, 

{  00  =  4>n+ i  =  o, 
nir /N,  C-  =  C. |_  and  C_  =  C+  =  0;  the  corresponding  frequency  is 


c ohh 


N 


20(2  —  cos(uh )  —  cos(o>/i)2) 
16  +  13  cos(c nh)  +  cos(tu/i)2 


(56) 


Figure  23  shows  the  discrete  frequencies  for  IV  =  9  degrees-of-freedom,  together 
with  the  graph  of  (56).  The  abscissa  is  the  scaled  mode  number  c oh/n  =  n/JV. 
Notice  that  the  present  scaling  of  the  mode  number  is  different  from  the  one  adopted 
for  linear  and  higher  order  p-methods  in  Sections  4.1.1  and  4.2. 


Remark  3  The  spurious  wave  components  e±luhA  do  not  contribute  to  the  normal 
modes  ( 55 ). 
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Fig.  23.  Analytically-computed  (discrete)  natural  frequencies  for  quadratic  /c-method 
(N  =  9). 

The  case  of  higher  order  elements  p  >  2  is  conceptually  similar,  though  the  higher 
is  p,  the  more  technical  the  analysis  becomes.  For  dispersion  analysis,  the  space  of 
discrete  solutions  (of  the  Helmholtz  equation)  has  dimension  2 p,  with  2p  —  2  spuri¬ 
ous  (linearly  independent)  solutions.  For  the  purpose  of  the  spectrum  analysis,  one 
has  to  split  the  eigenvalue/eigenvector  problem  into  boundary  and  interior  (stencil) 


2  We  emphasize  that  the  0^’s  are  the  control  variables  and  are  not  interpolated  by  the 
solution. 
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equations:  the  former  are  used  to  determine  which  solutions  of  the  interior  equa¬ 
tions  are  compatible  with  the  boundary  conditions,  thus  giving  admissible  normal 
modes.  They  have  a  structure  that  is  similar  to  (55)  but  the  last  p  —  1,  for  p  odd,  or 
p  —  2,  for  p  even,  correspond  to  evanescent  waves,  and  are  associated  to  “outlier 
frequencies”  (see  [7]). 

The  outlier  frequencies  disappear,  for  any  p,  when  a  nonlinear  parameterization  of 
the  domain  is  adopted,  through  the  uniform  distribution  of  control  points  shown 
in  Figure  4.  This  is  observed  numerically  but  a  sound  mathematical  explanation 
is  still  missing.  Indeed,  in  spectrum  analysis  the  /c-method  with  nonlinear  param¬ 
eterization  exhibits  a  more  complicated  behavior  than  the  /c-method  with  linear 
parameterization.  In  Figure  24  the  numerically  computed  ujh/co  (see  Remark  1) 
are  plotted  for  N  =  3,  10,  30  degrees-of-freedom:  it  is  possible  to  notice  that  the 
points  do  not  lie  on  an  underlying  curve  independent  of  N.  However,  they  converge 
towards  the  graph  of  the  analytical  relation  (56),  when  IV  — >  oo.  The  same  behavior 
is  observed  for  p  >  2  as  well;  the  computed  uh/u>  versus  the  scaled  mode  number 
converges,  for  N  — >  oo,  towards  the  analytical  relation  obtained  considering  only 
the  non-spurious  solutions  of  the  internal  stencil  equation. 


Fig.  24.  Numerical  spectrum  analysis  compared  to  the  analytically-computed  relation  (56) 
for  the  quadratic  /c-method  with  nonlinear  parameterization. 


5  Analytical  study  in  two  dimensions 


Spectrum  and  dispersion  analysis  can  be  put  in  relation  in  two,  and  more,  dimen¬ 
sions  as  well.  We  discuss  here  the  equivalence  between  the  study  of  propagation  of 
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numerical  waves  on  the  infinite  grid  and  the  study  of  discrete  natural  frequencies 
on  a  rectangular  domain.  At  the  discrete  level,  we  only  discuss  the  case  of  bilinear 
approximation,  but  the  same  concepts  extend  to  the  higher  order  p-method  (using 
the  techniques  of,  e.g.,  [9]),  and  fc-method  (dealing  with  the  boundary  equations  as 
in  the  one-dimensional  setting). 


Discretizing  the  Helmholtz  equation  by  bilinear  approximation  on  the  infinite  uni¬ 
form  mesh  of  square  elements  of  side  length  h,  the  discrete  equations  are 


—  811^^2)  +  m(Ai-1,A2-1)  +  m(Ai-1,A2)  +  m(Ai-1,A2+1) 

+  U(AltA2- 1)  +  U(A1,A2+ 1)  +  W(Ai+1,A2-1)  +  «(Ai+1,A2)  +  W(Ai+1,A2+1)) 

kh  1  1 

+  (4W(Ai,A2)  +  ^W(Ai-i,A2-1)  +  m(Ai-1,A2)  +  -M(Ai-1,A2+1)  +  m(Ai,A2-1) 

+  ti(A1,A2+l)  +  ^u(Ai+1,A2-1)  +  ^(A1+1,A2)  +  ^W(Ai+1,A2+1)),  VA  G  Z? , 

(57) 


where  uA  =  U(Ai,a2)  are  the  degrees-of-freedom  with  respect  to  the  standard  nodal 
basis. 


The  solutions  of  (57)  are  linear  combinations  of  discrete  plane  waves  uA  =  elhkh  A, 
where  (see,  e.g.,  [9])  the  discrete  wave-number  kh  =  ( k k%)  satisfies 


( khf 


—  ^  +  ^(/c/r)2^  (cos (k^h)  +  cos^/i)) 
)  cos (k^h)  cos(/c2  h)  =  0. 


(58) 


Likewise,  exact  plane-wave  solutions  of  (6)  are  w(x)  =  elk'x,  with  |k|  =  k. 


The  joint  plot  of  k hh  and  kh,  for  kh  =  1,  2,  3,  is  presented  in  Figure  25;  it  is  seen 
that  the  dispersion  error  depends  on  the  direction  of  propagation  of  the  wave.  A 
discrete  wave-number  kh  has  to  be  compared  with  the  corresponding  exact  wave- 
number  which  is  aligned,  that  is,  k  =  k-M^-r.  the  dispersion  error  is 


1  - 


W\ 


(F*IT 

.  The  plot  of  Tjpqj-  versus  khh  is  given  in  Figure  26. 


||k  — k|| 
||k^|| 


Turning  to  spectrum  analysis,  we  now  consider,  on  the  domain  =  (0,  L)  x  (0,1), 
a  mesh  made  of  nei  elements  per  direction.  Taking  into  account  the  boundary  con¬ 
ditions,  there  are  N  =  {nel  —  l)2  degrees-of-freedom.  The  discrete  equations  are 
analogues  of  (57),  by  invoking  the  change  of  notation  k  Uh  and  uA  =  (f>A.  The 
normal  modes  are 


(j)A  —  Csin^ihAi)  sm(uj2hA2),  VAi,  A2  =  0, . . . ,  neh  (59) 

with  1  <  ujiL/'k,  ui2L/'K  <  nei  —  1.  These  are  the  interpolants  of  the  exact  normal 
modes  0X  =  C  sin(cuiXi)  sin^a^)-  It  can  be  noticed  that  (59)  can  be  expanded  as 
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Fig.  25.  Numerical  wave-number  k''  and  exact  wave-number  k  for  k  =  1, 2,  3  in  two-di¬ 
mensions,  for  bilinear  element  approximation. 


Fig.  26.  Unified  dispersion  and  spectrum  analysis  in  two-dimensions,  for  bilinear  approxi¬ 
mation. 
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a  sum  of  four  plane  waves 


_  (J  fgi/i(wi,i^2)-A  _|_  gih(-o;i,-cJ2)'A  _  gih(— u;i,u;2)-A 


_iWu;i ojoVA  l 


(60) 


for  A  =  (Ai,  A2),  and  A1,  A2  =  0, . . . ,  nei.  Therefore,  we  are  exactly  in  the  sit¬ 
uation  of  the  previous  dispersion  analysis.  Fixing  toi  and  to2  by  changing  notation 
k\  uj\,  k2  (jJ-2,  and  using  (58),  the  discrete  frequency  is  obtained  from  the 
relation 


(J  -  \^h}l^2')  ~  (J  +  \^jhh)2)  (c°s(o;i/i)  +  cos(c o2h)) 

(A  (a jhh)2'  (61) 


+ 


cos(a;i  h)  cos(uo2h)  =  0. 


The  corresponding  exact  frequency  is  to  =  y to\  +  c o2.  Unifying  dispersion  and 

spectrum  analysis,  the  discrete-to-exact  frequency  ratio  —  versus  h,  u )2h  is  shown 
in  Figure  26. 

When  spectrum  analysis  is  carried  out  numerically,  as  in  [7]  for  example,  the  com¬ 
puted  frequencies  u )h  and  the  exact  ones  to  are  sorted  independently,  by  increasing 
magnitude,  and  then  associated.  In  one-dimension  this  produces  the  correct  associ¬ 
ation  of  vibrating  modes,  and  results  in  the  same  LOh /to  plot  obtained  analytically 
(see  Remark  1).  In  multi-dimensions,  however,  this  numerical  procedure  does  not 
guarantee  the  correct  association  of  modes  and  indeed  the  toh/to  numerically  cal¬ 
culated  differs  from  the  analytical  one.  However,  it  is  shown  in  Figure  27  that  the 
two  results  are  qualitatively  similar,  the  numerical  plot  presenting  less  oscillations 
but  still  revealing  the  correct  order  of  magnitude  of  the  error. 


6  Numerical  results 


In  this  section,  we  present  several  numerical  experiments  supporting  the  analytical 
results  previously  obtained.  In  particular,  for  both  k-  and  p-methods,  we  exam¬ 
ine  the  approximation  of  frequencies  and  modes  for  structural  dynamics  problems, 
while,  for  wave  propagation,  we  study  the  approximation  of  response  spectra,  as 
well  as  determine  solutions  of  a  ID  boundary  value  problem  proposed  in  the  liter¬ 
ature  (cf.  [22]). 

We  finally  show  preliminary  results  on  the  effects  of  under  integrating  NURBS 
discretizations,  compared  with  analogous  FEM  results. 
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Fig.  27.  Plots  of  ^  obtained  analytically  (by  making  first  the  correct  association  of  modes, 
and  then  ordering  them  such  that  ca  is  increasing)  and  numerically  (by  first  ordering  inde¬ 
pendently  the  frequencies  and  then  associating  them);  N  =  50. 


6. 1  Structural  vibrations 


We  present  results  of  frequency  calculations  for  ID  and  2D  cases,  and  we  study  the 
approximation  of  eigenmodes  in  ID. 


6.1.1  ID  probl  ems 

We  present  the  results  of  some  numerical  tests  performed  using  the  A  -mcthod  (i.e., 
NURBS)  and  p-method  (i.e.,  FEM)  on  a  ID  vibration  problem  (i.e.,  the  problem  of 
the  longitudinal  vibrations  of  an  elastic  rod).  In  the  previous  section,  this  has  been 
pointed  out  to  be  equivalent  to  the  classical  time-harmonic  dispersion  analysis  for 
ID  wave  propagation,  by  the  duality  principle. 

Before  comparing  NURBS  and  FEM  results,  we  briefly  remark  that  the  following 
plots  for  NURBS  are  obtained  using  a  nonlinear  parameterization  (as  described  in 
Section  3.3)  in  order  to  avoid  so-called  “outlier  frequencies”.  These  are  spurious 
frequencies  (or  discrete  optical  branches)  that  show  up  when  a  linear  parameteriza¬ 
tion  is  employed.  For  a  more  detailed  discussion  on  the  appearance  of  outliers  and 
how  to  eliminate  them,  the  reader  is  referred  to  Section  4.3  and  to  [7]. 

Figure  28  shows  a  comparison  of  k-  and  p-method  numerical  spectra  for  p  =  1, ...,  4 
(we  recall  that  for  p  =  1  the  two  methods  coincide).  Here,  the  superiority  of  the 
isogeometric  approach  is  evident,  as  one  can  see  that  optical  branches  of  spectra 
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diverge  with  p  for  classical  C°  finite  elements.  This  negative  result  shows  that  even 
higher-order  finite  elements  have  no  approximability  for  higher  modes  in  vibration 
analysis,  and  possibly  explains  the  fragility  of  higher-order  finite  element  methods 
in  nonlinear  and  dynamic  applications,  in  which  higher  modes  necessarily  partici¬ 
pate. 


n/N 


Fig.  28.  Comparison  of  A-mcthod  and  /> method  numerical  spectra. 

Finally,  following  [22],  we  study  the  approximation  of  eigenmodes  by  k-  and  p- 
methods.  In  [22],  Thompson  and  Pinsky  consider  the  ID  eigenproblem  correspond¬ 
ing  to  the  vibration  of  a  fixed-fixed  rod  of  unit  length,  discretized  with  21  degrees- 
of-freedom  (19  after  imposing  the  boundary  conditions)  and  with  quadratic  interpo¬ 
lations;  in  particular,  they  study  the  finite  element  approximation  of  the  18th  mode. 
In  order  to  compare  the  modal  approximation  properties  of  k-  and  p-methods,  we 
compute  eigenmodes  6,  9,  12,  15  and  18  for  this  same  problem  and  we  compare  the 
numerical  modes  with  the  analytical  ones,  namely, 

(p]{x)  =  sm(jnx),  (62) 

where  j  is  the  mode  number.  Figures  29-33  present  the  comparisons  and  clearly 
demonstrate  the  better  performance  of  the  ^-method  in  approximating  eigenmodes, 
especially  ones  corresponding  to  higher  frequencies. 


6.1.2  2D  problems 

We  conclude  with  a  comparison  of  k-  and  p-method  numerical  spectra  for  a  2D 
problem  (i.e.,  the  problem  of  the  transverse  vibration  of  a  membrane). 
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Fig.  29.  Exact  (red  dotted  line)  versus  numerical  (blue  solid  line)  6th  eigenmode.  Left: 
A:-mcthod  approximation;  right:  p-method  approximation. 


Fig.  30.  Exact  (red  dotted  line)  versus  numerical  (blue  solid  line)  9th  eigenmode.  Left: 
A:-method  approximation;  right:  p-method  approximation. 


Fig.  31.  Exact  (red  dotted  line)  versus  numerical  (blue  solid  line)  12th  eigenmode.  Left: 
A:-method  approximation;  right:  p-method  approximation. 


As  discussed  in  Section  5,  we  follow  [7]  and  represent  2D  spectra  obtained  numer¬ 
ically  as  we  did  for  ID,  that  is,  the  abscissae  are  the  normalized  numbers  of  modes 
sorted  from  the  smallest  to  the  highest  frequencies.  Figure  34  reports  the  numerical 
spectra  obtained  using  70x70  degrees-of-freedom.  The  results  exhibit  similarities 
to  the  ID  case  and  the  superiority  of  the  isogeometric  approach  is  also  clear.  Again, 
for  higher  frequencies,  finite  element  spectra  seem  to  diverge  with  p. 
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Fig.  32.  Exact  (red  dotted  line)  versus  numerical  (blue  solid  line)  15th  eigenmode.  Left: 
^-method  approximation;  right:  p-method  approximation. 


Fig.  33.  Exact  (red  dotted  line)  versus  numerical  (blue  solid  line)  18th  eigenmode.  Left: 
A:-method  approximation;  right:  p-method  approximation. 

6.2  Wave  propagation 


The  aim  of  this  section  is  to  compare  NURBS  elements  and  classical  finite  elements 
on  wave  propagation  problems.  In  particular,  we  study  the  problem  of  an  elastic 
rod  originally  proposed  in  [22]  and  we  use  the  k-  and  the  p-methods  to  compute 
the  numerical  frequency  response  spectra  and  to  solve  a  boundary  value  problem. 
We  note  in  passing  that  Figure  28  can  also  be  interpreted  as  representing  dispersion 
error  in  wave  propagation  according  to  the  duality  principle.  As  before  we  need 
to  make  the  interchange  ujh/cu  <-►  k / kh  and  khh  uh.  In  the  following,  all  the 
numerical  tests  are  carried  out  using  quadratic  and  cubic  elements. 


6.2.1  Frequency  response  spectra 

Following  [22],  we  start  from  the  governing  equation  of  the  steady-state  vibration 
problem  for  a  rod  of  length  L 


d  2(f> 
dx2 


F  kcj) 


0, 


(63) 
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Fig.  34.  Comparison  of  2D  A-mcthod  and  p-method  numerical  spectra.  Top:  entire  spec¬ 
trum.  Bottom:  detail  of  the  first  half  of  the  spectrum. 

with  boundary  conditions 

0(0)  =  0,  0(L)  =  0.  (64) 


The  solution  to  problem  (63)-(64)  can  be  written  as 


0(x,  k)  =  0 


sin  {k(L  —  x )) 
sin  (kL) 


(65) 


We  denote  by  0h(x.  k )  the  numerical  solution  for  the  discrete  methods.  Now,  the 
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dispersive  and  attenuation  characteristics  can  be  investigated  using  the  frequency 
response  function  for  both  k-  and  p-methods,  that  is,  we  compare  for  each  discrete 
method  the  values  of  \og10(R(x,  k))  with  the  corresponding  exact  values,  where 

R(x,  kh)  —  (j)h(x,  k)/(j).  (66) 


Figures  35  and  36  show  the  response  spectra  obtained  for  p  =  2  and  p  =  3  at 
x  =  L/ 10,  L/ 2,  and  91/ 10.  In  all  cases,  the  better  approximation  properties  of  the 
/c- method  are  evident,  as  well  as  the  very  poor  performance  of  the  p-method  within 
stopping  bands  (see  Figure  18  and  [22]  regarding  cubic  p-method  stopping  bands). 


6.2.2  Boundary  value  problem 

We  solve  the  ID  boundary  value  problem  for  different  choices  of  the  wave-number 
k  (taking,  e.g.,  </>  =  1  and  L  =  1).  In  order  to  have  meshes  with  elements  of 
the  same  length  (h  =  1/10)  independent  of  the  approximation  order,  we  use  21 
degrees-of-freedom  for  quadratics  and  31  for  cubics. 

In  Figures  37  and  38,  we  present  the  boundary  value  problem  results  for  both  k-  and 
p-methods  solved  with  k  =  10,  20,  30,  and  33  (i.e.,  within  the  p-method  stopping 
band)  for  quadratic  approximations,  and  with  k  =  10,  20,  30,  31.5  (i.e.,  within  the 
1st  p-method  stopping  band),  40,  50,  and  71  (i.e.,  within  the  2nd  stopping  band,  see 
Figure  18)  for  cubic  approximations.  In  the  case  of  NURBS,  no  evident  attenuation 
is  observed  within  the  1st  stopping  band,  which  is  very  narrow  and  has  a  very  small 
imaginary  part  (see  Figure  18).  The  phase  opposition  observed  in  the  p-method  for 
k  =  31.5  is  indeed  due  to  the  fact  that  at  k  =  lCU  <  31.5  an  exact  resonance 
peak  occurs,  which  is  approximated  by  the  p-method  slightly  after  k  =  31.5.  For 
the  sake  of  completeness,  Figure  39  shows  in  more  detail  what  happens  around  an 
exact  resonance.  We  wish  to  emphasize  that  these  resonance  peaks  do  not  appear 
in  the  frequency  response  spectra  of  Figure  36,  since  the  corresponding  (exact  and 
discrete)  eigenmodes  vanish  at  x  =  L/10,  L/2,  and  9L/10. 

These  figures  confirm  the  superiority  of  the  /.■-method  in  wave  propagation.  It  is 
noted  that  the  p-method  stopping  bands  result  in  spurious  attenuation  of  waves,  due 
to  imaginary  parts  of  discrete  wave-numbers. 


6.3  Under  integration 


We  perform  an  initiatory  investigation  of  approximate  integration,  an  issue  which 
is  of  considerable  importance  in  practical  analysis. 

For  p-  method  finite  elements,  numerical  quadrature  seems  to  be  fairly  well  under¬ 
stood.  If  we  assume  the  simple  case  of  integrating  polynomials  over  elements,  the 
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Fig.  35.  Frequency  response  spectra  for  p  =  2  computed  at  L/ 10  (top),  L/ 2  (middle),  and 
9L/10  (bottom). 

Gauss  rules  are  optimal  in  one  dimension,  and  often  utilized  for  tensor-product- 
based  multidimensional  elements.  For  NURBS  and  B-splines,  the  situation  does 
not  seem  to  be  very  well  understood,  even  in  one  dimension.  The  problem  here 
is  that  reduced  continuity  exists  at  all  knots.  Even  in  the  case  of  the  /c-method,  the 
continuity  across  knots  internal  to  a  patch  is  Ck  =  Cp~ 1 .  To  form  stiffness  matrices, 
we  need  to  differentiate,  and  form  products  of  derivatives,  resulting  in  polynomials 
of  order  2 (jp  —  1)  and  continuity  Up~2.  It  seems  what  is  required  for  B-splines  and 
NURBS  are  rules  that  account  for  the  degree  of  smoothness  across  knots,  and  the 
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Fig.  36.  Frequency  response  spectra  for  p  =  3  computed  at  L/10  (top),  L/2  (middle),  and 
973/10  (bottom). 

precise  basis  on  each  patch.  An  investigation  into  this  is  beyond  the  scope  of  this 
paper.  We  shall  use  an  approach  here  that  is  simple  and  effective,  but  very  ineffi¬ 
cient,  at  least  for  the  fc-method.  It  makes  use  of  the  observation  that  between  knots, 
NURBS  and  B-splines  are  C°°  and  so  are  their  derivatives.  Consequently,  Gauss 
rules  are  effective.  However,  this  amounts  to  overkill  compared  with  the  usual  C° 
finite  elements  because  their  basis  functions  are  C°°  across  internal  nodes  (i.e., 
knots).  Thus,  for  example,  in  one  dimension,  for  equal  order  k-  and  p-methods  one 
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k=10 


k=20 


Fig.  37.  Boundary  value  problem  solution  for  p  =  2  computed  with  k  =  10  (top-left), 
k  =  20  (top-right),  k  =  30  (bottom-left),  and  k  =  33  (bottom-right,  within  the  p-method 
stopping  band). 

would  be  using  p  times  as  many  points  for  the  /{-method  as  the  p-method,  because 
for  the  F- method  the  rule  needs  to  be  used  in  each  knot  interval,  whereas  for  the 
p-method  it  only  needs  to  be  used  for  each  element,  consisting  of  p  consecutive 
knot  intervals.  Nevertheless,  this  will  be  the  basis  of  this  initiatory  comparison.  It 
should  be  emphasized  that  this  means  that  many  more  points  are  being  used  in  the 
/{-method  than  for  the  p-method.  However,  conclusions  drawn  should  be  viewed  as 
preliminary,  at  least,  until  optimal  rules  are  developed  for  NURBS  and  B-splines. 


6.3.1  ID  spectrum  approximation 

We  start  by  considering  ID  problems  and  pointing  out  that  p  +  1  Gauss  points  are 
needed  in  order  to  exactly  integrate  both  mass  and  stiffness  matrices  obtained  from 
degree  p  basis  functions  (in  the  following,  we  will  refer  to  this  case  as  “full  inte¬ 
gration”).  Instead,  using  p  Gauss  points  (i.e.,  “under  integrating”  using  one  fewer 
Gauss  point),  the  mass  matrix  is  under  integrated  while  the  stiffness  is  still  ex¬ 
actly  integrated.  Using  less  than  p  Gauss  points,  we  under  integrate  both  mass  and 
stiffness.  We  remark  that  all  the  results  presented  prior  to  this  section  have  been 
obtained  using  full  integration. 
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k=20 


k=30 


k=31 .5 


k=40 


Fig.  38.  Boundary  value  problem  solution  for  p  =  3  computed  with  k  =  20  (top-left), 
k  =  30  (top-right),  k  =  31.5  (middle-left,  within  the  1st  p-method  stopping  band),  k  =  40 
(middle-right),  k  =  50  (bottom-left),  and  k  =  71  (bottom-right,  within  the  2nd  stopping 
band). 


We  first  study  what  happens  when  under  integrating  p-method  matrices  in  spectrum 
analysis.  Indeed,  in  this  case,  we  can  only  under  integrate  by  1  Gauss  point,  oth¬ 
erwise  stability  is  lost  (i.e.,  the  stiffness  matrix  becomes  singular).  Moreover,  the 
under  integrated  results  are  even  worse  than  the  fully  integrated  ones  in  that,  for 
fixed  p,  the  highest  frequency  error  diverges  as  the  mesh  is  refined.  See  Figures  40 
and  41  in  which  1000  control  points  were  used. 
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Fig.  39.  Boundary  value  problem  solution  for  p  =  3  computed  with  k  =  31.3  (top-left), 
k  =  31.4  (top-right),  k  =  31.5  (middle-left),  k  =  31.6  (middle-right),  k  =  31.7  (bot¬ 
tom-left),  and  k  =  31.8  (bottom-right),  illustrating  what  happens  around  the  exact  reso¬ 
nance  peak  occurring  at  k  =  107T. 

Better  results  are  obtained  under  integrating  ^-method  matrices  by  1  Gauss  point, 
as  shown  in  Figure  42.  Moreover,  it  is  interesting  to  observe  that  acceptable  results 
are  often  obtained  under  integrating  fc-method  matrices  by  even  more  than  1  Gauss 
points,  as  shown  in  Figure  43.  For  p  >  2  stability  is  always  lost  when  using  just  1 
Gauss  point,  so  in  the  tests  we  integrated  with  a  minimum  of  2  Gauss  points. 

Figure  44  shows  the  number  of  Gauss  points  needed  for  full  quadrature,  and  the 
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Fig.  40.  ID  numerical  spectra  for  linear  basis  functions  obtained  with  full  integration  com¬ 
pared  with  spectra  under  integrated  by  1  Gauss  point  (plotted  at  two  different  scales). 


Fig.  41.  ID  p-method  numerical  spectra  obtained  with  full  integration  compared  with  spec¬ 
tra  under  integrated  by  1  Gauss  point  (plotted  at  two  different  scales). 


Fig.  42.  ID  /c-method  numerical  spectra  obtained  with  full  integration  compared  with  spec¬ 
tra  under  integrated  by  1  Gauss  point  (plotted  at  two  different  scales). 

minimum  necessary  for  stability.  Perhaps  the  most  interesting  information  pre¬ 
sented  in  Figure  44  is  the  minimum  number  of  Gauss  points  needed  to  get  “accept¬ 
able”  results.  We  remark  that  the  “acceptable”  level  is  here  defined  by  a  subjective 
evaluation  of  the  spectrum  approximation  properties  given  on  the  basis  of  the  re¬ 
sults  reported  in  Figure  43  and  some  other  numerical  calculations.  It  is  interesting 
to  note  that  the  number  of  Gauss  points  needed  to  reach  acceptable  results  for  the 
/c-mcthod  is  described  by  the  expression  roundjup(p/ 2)  +  1  (where  round_up(-) 
is  the  round-toward-infinity  function).  Asymptotically,  the  slope  of  this  function  is 
1/2,  half  that  for  full  quadrature. 

A  rigorous  mathematical  explanation  of  the  effects  of  under  integration  on  p-  and 
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Fig.  43.  ID  /c-method  numerical  spectra  obtained  with  full  integration  compared  with  under 
integrated  ones,  for  different  choices  of  the  order  p. 
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Fig.  44.  Quadrature  limits  for  p-  and  /^-methods. 
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/c-mcthods  is  an  open  question.  Furthermore,  it  is  important  to  develop  new  and 
more  efficient  quadrature  rules  for  the  /c- method  so  that  its  full  potential  can  be 
reached.  The  stability  of  the  method  with  reduced  quadrature  suggests  to  us  that 
this  should  be  possible. 


7  Conclusions 


We  compared  the  approximation  properties  of  standard  C°  continuous  finite  ele¬ 
ments  with  NURBS  on  problems  of  structural  vibrations  and  wave  propagation. 
The  basis  of  the  comparison  is  the  same  number  of  degrees  of  freedom,  equiva¬ 
lently,  the  bandwith  of  the  matrix  system.  We  found  that  the  higher-modes  of  clas¬ 
sical  /(-method  finite  elements,  represented  by  so-called  “optical  branches”  of  the 
frequency  spectrum,  have  no  approximability  whatsoever  (not  surprising)  and  that 
the  errors  in  frequency  diverge  with  p  (very  surprising).  The  behavior  of  NURBS  is 
much  better.  The  entire  spectrum  converges  with/?.  This  suggests  to  us  that  NURBS 
present  the  possibility  of  higher  order  accuracy  and  robustness.  Heretofore,  within 
the  finite  element  method,  these  attributes  have  been  mutually  exclusive.  We  also 
articulated  a  “duality  principle”  which  provides  precise  correspondence  between 
spectrum  analysis  in  structural  dynamics  and  dispersion  analysis  in  wave  propa¬ 
gation.  Lastly,  we  performed  an  initial  study  of  reduced  quadrature,  being  fully 
aware  that  optimal  quadrature  rules  are  not  yet  available  for  NURBS.  Nevertheless, 
the  results  suggest  that  reducing  the  number  of  quadrature  points  for  NURBS  by  a 
significant  amount  is  feasible.  We  hope  to  pursue  this  issue  in  future  works. 
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